HISTÓRIA
Stanisław Marcin Ulam (1909-1984)
O matemático e físico polonês Stanisław Ulam, em uma tarde de 1946, entediado jogando Paciência, levantou a questão de quais seriam as probabilidades de um determinado jogador ganhar a partida. Então se dispôs a realizar vários cálculos, porém, percebeu que apenas a análise combinatória não seria o suficiente para chegar na resposta desejada.
Uma opção seria jogar várias vezes e ver quantas ganhou, entretanto, para atingir um número de jogadas suficientes ao ponto de obter a resposta levaria muito tempo, então teve a ideia de criar um algoritmo que simula quantas jogadas desejar. A computação da época era rudimentar, mas era o que tinha à disposição.
Por também ser físico e estar participando do Projeto Manhattan, percebeu que poderia aplicar sua ideia na Fissão Nuclear e explorar o comportamento das reações em cadeia de nêutrons. Assim fez e publicou o artigo “The Monte Carlo Method” em 1949.
Artigo “The Monte Carlo Method” - 1949
A SIMULAÇÃO
A Simulação de Monte Carlo, também chamada de Método de Monte Carlo, é um modelo probabilístico que gera amostras aleatórias prevendo resultados futuros com uma base escolhida. Em outras palavras, é uma simulação de cenários possíveis seguindo uma determinada função de densidade de probabilidade.
Weibull Inversa
A função de probabilidade escolhida para realizar a Simulação de Monte Carlo será a Weibull Inversa, mas também conhecida como Distribuição Fréchet, proposta por Jiang, Murthy e Ji (2001). O motivo da escolha foi sua massiva utilização em diversas áreas, principalmente na Medicina e Ecologia. A sua função densidade de probabilidade é esta a seguir:
\[ f(x) = \alpha\theta x^{-\theta-1}e^{-\alpha x^{-\theta}}, x > 0 \]
Sua implementação na linguagem R fica assim:
> dwi <- function(x, alpha, theta){
+ alpha*theta*x^(-theta-1)*exp(-alpha*x^(-theta))
+ }Sua função acumulada e sua implementação no R estão a seguir, respectivamente:
\[
F(x) = e^{-\alpha x^{-\theta}}, x > 0
\]
> pwi <- function(q, alpha, theta){
+ exp(-alpha*q^(-theta))
+ }Para a Simulação de Monte Carlo é necessário também obter a função quantil, pois, na simulação precisaremos gerar números aleatórios seguindo o nosso modelo probabilístico escolhido. Para encontrar a função inversa, precisaríamos considerar f(x) = q e isolar o x. Para alguns casos essa etapa seria um problema porque poderia acontecer de não existir uma solução analítica, todavia, para a Weibull Inversa não foi. Pensando em otimizar o processo, foi utilizado o software matemático Maplesoft.
O resultado obtido está a seguir, tanto a função quanto ela implementada no R:
> qwi <- function(p,alpha,theta)
+ {
+ exp(-log(-log(p)/alpha)/theta)
+ }\[ g(q)= e^{-log(-log(q)/alpha)/\theta}, x>0 \]
Para a simulação, precisaremos também de uma função que irá gerar números aleatórios segundo a distribuição escolhida. A chamaremos de ‘rwi’.
> rwi <- function(n, alpha, theta)
+ {
+ u <- runif(n)
+ rv <- qwi(u, alpha, theta)
+ return(rv)
+ }Existem ainda duas funções que são necessários implementar antes de iniciarmos de fato a simulação: a Função Log Verossimilhança e o Estimador de Máxima Verossimilhança. Em relação ao primeiro item, basta fazermos o somatório do log da Função Densidade de Probabilidade. Segue sua implementação:
\[ l(\alpha,\theta;x) = \sum\limits_{i=1}^{\mbox{n}}log(\alpha\theta xi^{-\theta-1}e^{-\alpha xi^{-\theta}}) \]
> ll <- function(x, par)
+ {
+ sum(log(dwi(x,alpha = par[1],theta = par[2])))
+
+ }A segunda função que está faltando é a do Estimador de Máxima Verossimilhança. Se fôssemos calcular na mão, precisaríamos utilizar a matriz Hessiana para encontrar, contudo, computacionalmente temos algumas opções, dentre elas o pacote ‘Fitdistrplus’ e o ‘Optim’. O ‘Fitdistrplus’ tem sua vantagem, pois, para determinados problemas costuma ser mais rápido, todavia, acabei optando por usar o ‘Optim’ por já fazer parte da base da Linguagem R. Sua implementação está a seguir:
> emv.weibull <- function(x, par){
+
+ fit <- optim(par = c(par[1], par[2]), fn = ll, x = x,
+ control = list(fnscale = -1))$par
+ if(!is.numeric(fit)) fit <- NA
+ return(fit)
+ }Agora irei criar uma semente para ser possível replicar os resultados que virão a seguir e também determinarei os parâmetros que irei simular. A escolha dos parâmetros foram baseados no artigo ‘Modifications of the Weibull distribution: A review’, material que consultei para definir as Funções Densidade e Acumulada da Weibull Inversa. É necessário também criar uma variável que contenha todos os cenários possíveis permutados entre os parâmetros escolhidos. Para a tarefa, usarei o ‘expand.grid’. A seguir estarão todos os cenários possíveis.
> # Estrutura da Simulação
>
> set.seed(16011995) # Definir a semente para resultados reprodutíveis
>
> alpha <- seq(1, 3, 1) # Definir os valores paramétricos
> theta <- seq(1, 3, 1)
>
> param <- expand.grid(alpha, theta) # CenáriosPossíveis Cenários
| Alpha | Theta | |
|---|---|---|
| 1 | 1 | 1 |
| 2 | 2 | 1 |
| 3 | 3 | 1 |
| 4 | 1 | 2 |
| 5 | 2 | 2 |
| 6 | 3 | 2 |
| 7 | 1 | 3 |
| 8 | 2 | 3 |
| 9 | 3 | 3 |
Agora iremos determinar a quantidade de simulações, os tamanhos de amostras possíveis e quatro matrizes que armazenarão os Viéses e os Erros Quadráticos Médios de cada parâmetro.
> B <- 10 # Número de simulações
>
> nmax <- 100 # Tamanho de amostra máximo
> enes <- seq(10, nmax, 10) # Tamanhos de amostras
>
> v.alpha <- matrix(nrow = length(enes), ncol = nrow(param))
> # Viés do estimador do parâmetro alpha
>
> v.theta <- v.alpha # Viés do estimador do parâmetro theta
>
> e.alpha <- matrix(nrow = length(enes), ncol = nrow(param))
> # Erro-quadrático-médio do estimador do parâmetro alpha
>
> e.theta <- e.alpha # Erro-quadrático-médio do estimador thetaAgora iniciaremos a simulação. O código implementado a seguir irá calcular o Estimador de Máxima Verossimilhança para cada x[i] de todos os tamanhos de amostras possíveis em cada cenário contido na variável ‘param’. Exemplo: o cenário 1, que possui \(\alpha = 1\) e \(\theta = 1\), a simulação irá calcular o EMV de todos os x[i]‘s das amostras de tamanho 10, 20, 30, 40, 50, 60, 70, 80, 90 e 100. Depois, armazenarão em ’v.alpha’, ‘v.theta’, ‘e.alpha’ e ‘e.theta’ os viéses e EQM’s dos parâmetros.
Dessa forma, conseguiremos definir qual é o conjunto de parâmetros que possui o menor Viés e o menor Erro Quadrático Médio.
> # Simulação
>
> for(i in 1:nrow(param)){
+ k <- 1
+ X <- rwi(nmax * B, alpha = param[i,1], theta = param[i,2])
+ X <- matrix(X, ncol = B, nrow = nmax)
+
+ for(n in enes)
+ {
+ x <- data.frame(X[1:n,])
+ fit <- sapply(x, emv.weibull, par = param[i,])
+
+ v.alpha[k,i] <- mean(fit[][1,] - param[i,1], na.rm = TRUE)
+ v.theta[k,i] <- mean(fit[][2,] - param[i,2], na.rm = TRUE)
+
+ e.alpha[k,i] <- mean((fit[][1,] - param[i,1])^2, na.rm = TRUE)
+ e.theta[k,i] <- mean((fit[][2,] - param[i,2])^2, na.rm = TRUE)
+
+ k <- k + 1
+
+ cat(i, n, param[i,1], param[i,2], "\n")
+ }
+ }Após iniciar a compilação do código implementado, esse processo pode levar de alguns segundos a algumas horas dependendo da quantidade de simulações que for realizar. No meu caso, se colocarmos para fazer 100 mil simulações, irá levar aproximadamente 6 horas para finalizar. O tempo longo é um ponto negativo, porém a seguir vamos entender o porquê de usar um número alto assim.
Mas antes de irmos para os resultados, após a compilação, o ideal é plotar os gráficos e armazená-los para não correr o risco de ter que rodar novamente as linhas de comando. Os comandos para realizar esta tarefa estão a seguir:
> png(filename = "viesalpha10.png",width = 25, height = 21,
+ units = "cm",res = 1200)
> matplot(v.alpha, type = 'b', xlab = 'Amostra', ylab = NULL,
+ main = 'Viés do Estimador do Parâmetro Alpha - 10 Simulações')
> dev.off()
>
>
> png(filename = "viestheta10.png",width = 25, height = 21,
+ units = "cm",res = 1200)
> matplot(v.theta, type = 'b', xlab = 'Amostra', ylab = NULL,
+ main = 'Viés do Estimador do Parâmetro Theta - 10 Simulações')
> dev.off()
>
>
> png(filename = "eqmalpha10.png",width = 25, height = 21,
+ units = "cm",res = 1200)
> matplot(e.alpha, type = 'b', xlab = 'Amostra', ylab = NULL,
+ main = 'EQM do Estimador do Parâmetro Alpha - 10 Simulações')
> dev.off()
>
>
> png(filename = "eqmtheta10.png",width = 25, height = 21,
+ units = "cm",res = 1200)
> matplot(e.theta, type = 'b', xlab = 'Amostra', ylab = NULL,
+ main = 'EQM do Estimador do Parâmetro Theta - 10 Simulações')
> dev.off()RESULTADOS
A seguir mostraremos os resultados das 10, 100, 1000, 10 mil e 100 mil simulações. Lembrando que no eixo das abscissas serão representadas as amostras numeradas de 1 até 10, porém significa que a amostra 1 tem 10 observações, a 2 tem 20 e assim por diante até chegar na amostra de tamanho 100, como foi definido no script.
10 simulações
100 simulações
1000 simulações
10 mil simulações
100 mil simulações
